function reed_analysis()
%
% analyze the reed data sets.

% handle multiple selections first
global CONTROL

sf = getmainselection;
frec = 100000;
lrec = 0;
if(sf > 0)
    pflag = getplotflag;
    QueMessage('Reed-Mitchell Analysis', 1); % clear the que
    h = findobj('Tag', 'Reed-Mitchell'); % check for pre-existing window
    if(isempty(h)) % if none, make one
        h = figure('Tag', 'Reed-Mitchell', ...
            'Units', 'normalized', ...
            'MenuBar', 'figure', ...
            'WindowButtonMotionFcn', 'datac(''mouse_motion'', gcbf);', ...
            'WindowButtonDownFcn', 'datac(''mouse_down'', gcbf);', ...
            'WindowButtonUpFcn', 'datac(''mouse_up'', gcbf);', ...
            'NumberTitle', 'off');
        datac('addwindow', 'Reed-Mitchell');

    end
    figure(h); % otherwise, select it
    clf; % always clear the window...
    %Command Menu
    % uimenu('Label', 'Close &Window', 'Position', 1, 'callback', 'close(findobj(''Tag'', ''CC_IV''));' );
    %uimenu('Label', '&Print', 'Callback', 'print;');
    % uimenu('Label', 'Print&Setup', 'Callback', 'printdlg;');
    tipi = [];
    pipi = [];
    gtseq = {};
    gpseq = {};
    for i = 1:length(sf)
        [t, p, fr, lr, tseq, pseq] = ra_2(sf(i), pflag);
        if(fr < frec)
            frec = fr;
        end;
        if(lr > lrec)
            lrec = lr;
        end;

        tipi = [tipi t];
        pipi = [pipi p];
        gtseq{i} = tseq;
        gpseq{i} = pseq;
        
    end;

    h = findobj('Tag', 'Reed-Mitchell-summary'); % check for pre-existing window
    if(isempty(h)) % if none, make one
        h = figure('Tag', 'Reed-Mitchell-summary', ...
            'Units', 'normalized', ...
            'MenuBar', 'figure', ...
            'WindowButtonMotionFcn', 'datac(''mouse_motion'', gcbf);', ...
            'WindowButtonDownFcn', 'datac(''mouse_down'', gcbf);', ...
            'WindowButtonUpFcn', 'datac(''mouse_up'', gcbf);', ...
            'NumberTitle', 'off');
        datac('addwindow', 'Reed-Mitchell-summary');

    end
    figure(h); % otherwise, select it
    clf; % always clear the window...
    indexpos = 4;
    sym = {'kv', 'go', 'bs', 'ro', 'yp', 'ko', 'gs', 'bd', 'rp', 'yv'};
    colrs = ['k', 'g', 'b', 'r', 'y', 'k', 'g', 'b', 'r', 'y'];
    set(h, 'Name', sprintf('Reed-Mitchell-Summary - File: %s', CONTROL(sf(1)).filename))
    figure(h);
    subplot('Position', [0.1, 0.55, 0.8, 0.4]);
    plot(tipi, pipi, [sym{indexpos} '-']);
    set(gca, 'YLim', [0 1.2]);
    set(gca, 'Xlim', [0 1.2*max(tipi)]);
    xlabel('Time pulses 1-3 (msec)');
    ylabel('Mean Firing Prob (0-1)');
    title(sprintf('R-M: File %s', CONTROL(sf(1)).filename));
    subplot('Position',[0.1,0.05,0.8,0.45])
    for i = 1:length(sf)
        sp = find(gpseq{i} == 1);
        nosp = find(gpseq{i} == 0);
%        plot(gtseq{i}(sp), ones(size(gpseq{i}(sp)))*i, ...
%            char(sym(mod(i,10)+1)), 'markerfacecolor', colrs(mod(i,10)+1));
       maxt = max([gtseq{i}(:)])
       plot(gtseq{i}(sp), ones(size(gpseq{i}(sp)))*tipi(i), ...
           sym{indexpos}, 'markerfacecolor', colrs(indexpos));
       hold on;
%        plot(gtseq{i}(nosp), ones(size(gpseq{i}(nosp)))*i, char(sym(mod(i,10)+1)));
        plot(gtseq{i}(nosp), ones(size(gpseq{i}(nosp)))*tipi(i), sym{indexpos});
    end;
    set(gca, 'Ylim', [0 1.2*max(tipi)]);
    set(gca, 'xlim', [0 1.2*max(tipi)]);
    xlabel('Position of 2^{nd} pulse, msec');
    ylabel('1^{st}-3^{rd} interval (msec)');
    
    %axis([0,1,0,1])
    %axis('off')
    %text(0,0.2,sprintf('%-12s R[%d:%d] %-8s',CONTROL(sf(1)).filename, frec, lrec, CONTROL(sf(1)).protocol), 'Fontsize', 8);
end;

function [tipi, pipi, fr, lr, tseq, pseq] = ra_2(sf, pflag)
% sf selects the record we will analyze this call
% pflag is the plot flag... as checked on the main display window.
global VOLTAGE DFILE CONTROL SFILE

tipi = [];
pipi = [];

[DFILE, err] = analysis_setup(DFILE, sf);

if(err ~= 0)
    return;
end;
fsize = 7;
msize = 3;
msg_pos = [0.37 0.00 0.15 0.07];

h = findobj('Tag', 'Reed-Mitchell'); % check for pre-existing window
if(isempty(h)) % if none, make one
    QueMessage('Can''t make Reed-Mitchell figure window', 1);
    return;
end;
figure(h);
set(h, 'Name', sprintf('Reed-Mitchell - File: %s', CONTROL(sf).filename))
figure(h);

tb=make_time(DFILE);
plot(tb', VOLTAGE');
hold on;

dt = DFILE.rate(1)*DFILE.nr_channel(1)/1000;
pars.rise=0.15;
pars.decay=0.5;
pars.threshold=3;
pars.sign=1;
pars.dispflag = 0;
pars.lpfilter=10000;
pars.template_type=1;
spikethr = -45;

x=spike_shape(pars); % get the spike times...

% now get the time window information for this block
[sfile, df] = block_info([CONTROL(sf).path CONTROL(sf).filename], DFILE.dblock);
tpw13 = number_arg(sfile.IPI.v); % interval from first to last pulse
tpwr = seq_parse(sfile.Sequence.v); % get roving pulse intervals.
t0 = number_arg(sfile.Delay.v);
tseq = tpwr{:}-t0;
pseq = zeros(size(tseq));
figure(h);
nsp = 0;
nseq = length(tseq);
if(length(x) < nseq)
    nseq = length(x); % may be early termination
end;
for i = 1:nseq % for all of the records in the sequence...
    xn = [x{i}]; % convert from cell to simple array
    if(~isnan(xn)) % make sure there are spikes
        for k = 1:length(xn)
            xs = xn(k);
            xt=number_arg(dt*xs);
            c=curvature(tb(i,:), VOLTAGE(i,:), xs, 2);
            if((VOLTAGE(i, xs) >= spikethr) & (c < 2))
                plot(tb(i,xs), VOLTAGE(i,xs), 'go', 'markerfacecolor', 'g' );
                fprintf(1, '%2d.%d %.4f  :  %.4f  %.3f\n', i, k, tseq(i), xt, c);
                %            [tb(i,xs), VOLTAGE(i,xs), c]
                nsp = nsp + 1; % count up spikes
                pseq(i) = 1;
            else
                fprintf(1, '%2d.%d %.4f  %.3f  :  event at %.3f ms failed criteria\n', i, k, tseq(i), c, xt);
                plot(tb(i,xs), VOLTAGE(i,xs), 'rs', 'markerfacecolor', 'r' );
            end;
        end;
    else
        fprintf(1, '%2d. %.4f :   no spike\n', i, tseq(i));
    end;
end;
fprintf('Firing Probability for 1-3 IPI of %.4f ms is %5.3f\n', tpw13, 100*nsp/nseq);
tipi = tpw13;
pipi = nsp/nseq;
fr = DFILE.frec;
lr = DFILE.lrec;
return;


function [allev] = spike_shape(pars)
% get spikes from traces using template matching - this is used to get
% spikes from an artifact ridden trace.
% 8/2/04 P. Manis
%

global VOLTAGE DFILE

samplerate =  (DFILE.rate(1)*DFILE.nr_channel(1)/1000);
rise=pars.rise;
decay=pars.decay;
threshold=pars.threshold;
sign=pars.sign;
dispflag = pars.dispflag;
lpfilter=pars.lpfilter;
template_type=pars.template_type;

vmean=[];
allev=cell(1, size(VOLTAGE, 1));
for k = 1:size(VOLTAGE,1)
    vmean(k) = mean(mean(VOLTAGE(k,1:10)));
    [evn, isamp, icoff, crit, template, predelay] = ClementsBekkers(VOLTAGE(k,:)-vmean(k), samplerate, rise, decay, threshold, sign, dispflag, lpfilter, template_type);
    if(~isempty(evn))
        [tmax, imax] = max(template);
        allev{k} = evn+imax-floor(predelay/samplerate)+2; % store result
    end;
end;

return;


function [c]=curvature(x, y, k, m)
% calculate the curvature of the trace v at the point k
% return the result in c.
% use to find the tops of spikes and other kinds of events.
% m is the width (in data points) to take around the point k.
%

dy=diff(y);
d2y=diff(dy);
nu = abs(mean(d2y((k-m-2:k+m-2))));
de = (1+mean(dy(k-m-1:k+m-1))^2)^(3/2);
c = nu/de;
return;


